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Abstract 

We study the adaptive dynamics of predator-prey systems modeled by a dy- 
namical system in which the traits of predators and prey are allowed to evolve 
by small mutations. When only the prey are allowed to evolve, and the size of 
the mutational change tends to 0, the system does not exhibit long term prey 
coexistence and the trait of the resident prey type converges to the solution of 
an ODE. When only the predators are allowed to evolve, coexistence of preda- 
tors occurs. In this case, depending on the parameters being varied, we see (i) 
the number of coexisting predators remains tight and the differences in traits 
from a reference species converge in distribution to a limit, or (ii) the number 
of coexisting predators tends to infinity, and we calculate the asymptotic rate 
at which the traits of the least and most "fit" predators in the population 
increase. This last result is obtained by comparison with a branching random 
walk killed to the left of a linear boundary and a finite branching-selection 
particle system. 
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1. Introduction 

The rapidly developing field of adaptive dynamics emphasizes the com- 
bined effects of evolution and ecological interactions on population dynamics. 
To describe the general framework of this theory, consider a population of 
individuals, each attached with a trait or strategy x that characterizes its 
ability to survive and propagate. The current distribution of traits governs 
the population dynamics by describing interactions between different individ- 



* Corresponding author: Tel.: +1 607 255 8262; fax: +1 607 255 7149 
Email addresses: rtdl@cornell.edu (Rick Durrett), jm858@cornell.edu (John 
Mayberry) 

1 Partially supported by NSF grant DMS 0704996 from the probability program. 
2 Partially supported by NSF RTG grant DMS 0739164. 



Preprint submitted to Elsevier 



March 23, 2010 



uals and their environment. Underlying this ecological process, is a slower 
mutational process that occasionally introduces a new trait into the popula- 
tion. Understanding the interplay of these fast ecological and slow mutational 
processes is the primary objective of adaptive dynamics. Foundations of the 
theory were laid in early 1990's by Hofbauer and Sigmund [22], Metz, Nis- 



bet, and Geritz [26], and Dieckmann, Marrow, and Law [14J and focused on 



macroscopic models, i.e., ODE models describing large population limits. In 



their 1996 paper [13] . Dieckmann and Law suggest that: "A proper mathe- 
matical theory of evolution should be dynamical... The dynamics ought to be 
underpinned by a microscopic theory" . A rigorous foundation for microscopic 
models via multi-type branching processes has now been developed (see for ex- 
ample Champagnat and Lambert 1(| and Champagnat, Ferriere, and Meleard 



)• 

In this paper, we study the dynamics of coexistence that arise as a con- 
sequence of introducing rare, small mutations into a model for predator-prey 
interactions. The novelty of our work lies in the establishment of coexistence 
of a large number of types. This phenomenon is known as polymorphic evo- 
lution. Two other notable examples of polymorphic evolution in the adaptive 
dynamics literature are: (i) evolutionary branching (first described in Geritz 



et al. 21] and more recently studied in Champagnat and Meleard Uj from a 



microscopic perspective) which describes coexistence of types with diverging 



traits and (ii) the Tube Theorem of Geritz et al. 20], [19] where coexistence 
of a resident and invading type with similar survival strategies occurs inside of 
a "tube" in which the sum of the invader and resident population sizes stays 
close to the former resident attractor. The situation we encounter more closely 
resembles the second scenario: types with very similar traits can coexist. In 
our model, this is due to the fact that interspecific competition (competition 
with individuals of different types) has less of an effect than intra-specific 
competition (competition with individuals of the same type). 

Since our focus will be on the dynamics of the random process of types that 
emerges from our underlying mutational process, we shall take a macroscopic 
perspective of population dynamics, using a Lotka-Volterra system of ODE's 
to describe predator-prey interactions. In particular, we will suppose that if 
we have M prey types and N predator types, then the densities are governed 
by the ODE's 



dui 

-dt = 



^ = V 3 | a iJ U i ~ 6 J ~ V ) ( L1 ) 

where the Ui, 1 < % < M are the densities of the prey, and the Vj, 1 < 



2 



j < N are the densities of the predators. Our main interest is the effect 
of small mutations in the resident types on the equilibrium behavior of (11. ip . 
While the co-evolutionary case in which both predator and prey are allowed to 
vary is certainly of interest and can lead to exotic behavior (see, for example, 



Dieckmann et al. [14j and Dercole et al. |12j). we will here only consider the two 



cases of fixed predator /evolving prey and evolving predator/ fixed prey. Such 
examples are also of interest and have been studied in laboratory experiments 
(see, for example, Jones and Ellner 

Following the usual approach in adaptive dynamics, we shall assume that 
mutations take place on a much slower time scale than the population dynam- 
ics reach equilibrium. In particular, suppose that we are considering predator 
evolution and we currently have k predator types (which we call the resident 
types) and one prey type coexisting in equilibrium. We introduce a small 
density of a new type of predator, called the mutant type, with traits chosen 
according to a specified mutation distribution, and let the densities evolve ac- 
cording to (11.11) until a possibly new equilibrium is reached before introducing 
the next mutant type into the population. By traits, we mean the parameters 
in (II. ip (birth, death, and consumption rates) that characterize each individ- 
uals ability to survive and propagate. We then repeat this process, using only 
those predators which could coexist at the previous step. In this way, we 
obtain a Markov chain of resident types with transitions determined by the 
outcomes of the ecological interactions between the previous residents and the 
mutant. Once we have introduced some preliminary results for the ODE, we 
will formulate this process more precisely. 

Our evolutionary algorithm leads to a variety of different scenario depend- 
ing on the underlying parameters being varied. 

(a) If we have a single, fixed predator and allow mutations in the prey's traits 

(a, (3), then prolonged coexistence of prey does not occur and the traits 
of the resident prey converge to the solution of an ODE (see Theorem 

(b) If we have a single, fixed prey and allow the consumption rate a of preda- 

tors to evolve, then coexistence of predators occurs, but the number of 
coexisting predators remains tight and the differences of the parameters 
from a reference type converge in distribution to a limit (see Theorem 

[2D. 

(c) If we have a single, fixed prey and allow the death rate S of predators to 

evolve, the number of coexisting predators tends to infinity and we can 
calculate the speed at which the traits of the least and most fit predators 
in the population increase (see Theorem [3]). 
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In all three cases, our results are more mathematically interesting than biolog- 
ically relevant since in (b), for example, the consumption rate of all predators 
currently present in the population increases without bound. 

The remainder of this section is dedicated to statements of these results 
and some conjectures for future research. Proofs of the three results in (a) - 
(c) above will be contained in Sections [2J15J 

1.1. Prey Evolution 

We begin with the case in which we have a single predator with fixed death 
rate 5 > and we allow prey types to evolve. Throughout the remainder of 
this subsection when we refer to (11 .ip , we shall always assume that N — 1 and 
let v denote the density of our single predator. Prey types are characterized 
by their two dimensional trait vectors y = (a, /3) 6l 2 and we say that prey 
types yi, ...,yM can coexist with the predator if whenever we run ( II. ip started 
from positive initial densities, v(0),Ui(0) > 0, the densities v(t), Ui(t) remain 
bounded away from zero for all time. Our first step is to discuss criteria for 
coexistence. We use the notation u = (u\,U2, ...,Um) and 

r M ,i = {(u, v) e R M+1 -.v, Ul >o, < 1} 

r Afti = {( u ' v ) £ r M,i : v,Ui > 0, Mi E J} 

for any J C {1, 2, ...M}. If J = {1, M }, we simply write = T J ^ V Note 
that Tm,i is invariant under ( II. ip . Here and elsewhere, we shall use | • | to 
denote the cardinality of a finite set, the absolute value of a real number, and 
the Lebesgue measure of a set of real numbers depending on the context. 

Proposition 1.1. For all prey yi e R 2 , i < M, with different birth rates, ( 11.11) 
has an explicitly calculable equilibrium o = (ox, cr M , (Tm+i) £ ^m,i which is 
globally attracting on V^i- Furthermore, if J a = {i < M : <jj > 0} ; then 
\J<r\ < 2 and a is globally attracting on as well. 

Proposition 11.11 follows from Lemmas 12.11 and 12.21 in Section [2j Since we 
need explicit formulas for the equilibria of (ll.ip . we will need to redo some 
standard results on Lotka-Volterra equations (see, for example, Chapters 13 
and 15 of Hofbauer and Sigmund 23] and Chapter 3 in Takeuchi 28|) to prove 



our results. 

We are now ready to describe the Prey Evolutionary Process (Prey EP). 
This process is a continuous time Markov jump process which keeps track of 
the current resident prey types in the population. Proposition II. II tells us that 
we will never have more than two coexisting prey types so at time t, the state 
of the Prey EP is Y(t) = (Y 1 (t),Y 2 (t)) E R 2 x 1R 2 . For initial conditions, we 
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set ^2(0) = (0,0), i.e., the second prey species is initially absent, and choose 
any 1^(0) = (a(0),/3(0)) satisfying 



^(0) E V = Ua,(3) e R\ : /3 > — ^ > 1 J . 

The reason for this choice of Yi(0) is that if M = 1, the globally attracting 
equilibrium described in Proposition 11.11 satisfies 01,02 > if and only if the 
prey type has traits in V (see (12.31) in Section [2]). As long as Yi(t) ^ (0,0), 
mutational events occur at rate 1 and after a mutation, the transitions for 
Y(t) are determined by the following procedure. We pick one of the non- 
zero Yi(t—), i = 1,2 at random and choose Y new = (a new , (3 new ) uniformly 
from B e (Yi(t-)), the ball of radius e around !*(*-). If Y 2 (t-) ^ (0,0), let 
a be the equilibrium obtained in Proposition 11.11 when M = 3 and the prey 
have traits 3/1,3/2,3/3 = Y 1 (t-),Y 2 {t-),Y new . If Y 2 (t-) = (0,0), then let o be 
the equilibrium obtained in Proposition 11.11 when M = 2 and the prey have 
traits 3/1,3/2 = Yi(t—),Y new . Note that since the probability of inserting a 
mutant with the same birth rate as one of the residents is 0, we do not have 
to worry about the exceptional case where Proposition 11.11 does not apply. If 
I Jo- 1 = 2, then we set Yi(t) and ^2^) equal to the parameter values of the 
two prey with positive equilibrium densities. If \J a \ = 1, then we set Yi(t) 
equal to the parameter values of the single prey with positive equilibrium 
density and take Y 2 (t) = (0, 0). If | J a \ = 0, we set Yi(t), Y 2 (t) = (0, 0) and the 
process enters an absorbing state. We say that the population is monomorphic 
when Y 2 (t) = (0,0) and refer to the events where Y 2 (t) jumps from (0,0) as 
coexistence events. 

Our first main result says that in the small mutation limit, the population 
is essentially monomorphic. 

Theorem 1. Let T > 0. As £ -)■ ; Y((t) = Yi(t/e) ->■ yi(t) in probability 
uniformly on [0, T]. 3/1 (t) is the unique solution to the ODE 

with initial conditions 3/1(0) = Yi(0) where A/"(-) is explicitly calculable, see 
\2.1$ . Furthermore, if we let Y 2 £ (t) = Y 2 {t/e) and 

N t £ = \{s<t:Y 2 £ (s-)=0,Yi(s)^0}\ 

t < T, be the number of times Y 2 jumps from before time t/e, then as e — > 0, 

N £ => N 

where N is a nonhomogeneous Poisson Process on [0, T] and " =$■" denotes 
convergence in distribution. 
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We will prove Theorem [T] in Section [2j The proof reveals that \{t < T : 
Y 2 £ (£ At) > 0}| — > (see (I2.8P ) which justifies our earlier claim that prolonged 
coexistence of prey does not occur. The constant on the right hand side of 
( 11.21) is EY + when (X, Y) is chosen at random from the ball of radius 1 
and appears due to our choice of mutation distribution. (11.21) is essentially a 
special case of the "Canonical Equation of Adaptive Dynamics," see (6.2) in 



Dieckmann and Law [13| or (1) in Champagnat and Lambert [10(]. In words, 
evolution proceeds in the direction of fastest increase in fitness. We do not 
have an explicitly defined fitness, but the infinitesimal drift in the traits is 
perpendicular to the region of values that cannot invade the resident. 

The limiting ODE in Theorem [1] is not biologically sensible because the 
prey birth rate increases without bound. This could be remedied by restricting 
the permitted values of (a,/3) to a curve, or making them functions of other 
underlying parameters (see for example Dieckmann et al. [3], where traits 
depend on the "body size" of the predator and prey). However, since our 
main interest in including Theorem [1] is as a contrast to the results on predator 
evolution below, we do not here endeavor to carry out the details of this more 
realistic scenario. 



1.2. Predator Evolution 

We now consider the case where we have a single prey with fixed birth 
rate > 1 and density u (M — 1 in ( 11.11) ). but we allow our predators 
to evolve. Predators are characterized by their two dimensional trait vector 
x = (a, 5) G R+. As in the previous section, the first step is to develop 
a criteria for determining coexistence of multiple predators. The next result, 
which is proved in Section [31 tells us that this can be done by checking a simple 
algebraic condition. In order to state the result, we define the characteristic 
ratio of predator x = (a, 5) as £ = a/ 5 and use the notation 

r 1)JV = {(u, v h v N ) G R N+1 : < u < 1, Vi > 0} 
T+ N = {(u, v h v N ) G R^ +1 : < u < 1, Vi > 0}. 

Note that is invariant for (11.11) . 

Proposition 1.2. For any N > 1, ( 11.11) has a unique, globally attracting 
equilibrium a for initial densities in T^ N . This equilibrium has the following 
characterization. Suppose that predators X\,...,x^ are ordered by increasing 
characteristic ratios. Then the globally attracting equilibrium has a positive i th 
component if and only if i = l,...,m where m < N is the largest value of k 
satisfying the condition: 

k 

<r-Ph (1-3) 

i=i 
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and r = f3 — 1 > is the intrinsic growth rate of the prey. 



We prove Proposition 11.21 in Section [3j The definition for the Predator 
Evolution Process (Predator EP) is similar to the definition of the Prey EP 
except the state space is now (M 2 ) N as there is no limit on the number of 
predators that can coexist. The state of the process at time t is 

X(t) ^ ((a 1 (t),S 1 (t)),(a 2 {t),6 2 (t)),...). 

For initial conditions, we choose («i(0), #i(0)) so that a 1 (0)r/ J 5 > #i(0) and 
set afc(0),/3fe(0) = for all k > 2. Our choice of ai(0), <5i(0) guarantees 
that the initial predator can coexist with the prey (see Section |3]). As long 
as (ai(t) , j3x(t)) ^ (0,0), mutations occur at rate 1 and if a mutation oc- 
curs at time t, transitions are described by the following rules. Define N t = 
max{i : (aci(t— ), 0i(t— )) ^ (0,0)}, choose one of (aa(t— ), 6i(t— )), i < N t at 
random, and introduce a new mutant predator with traits a new = c*i + eU\ 
and 5 new = Sie eU2 where U\, U% ~ Uniform [-1,1]. We then order the predators 
(ai(t-),/3 2 (t-)), • • • (a Nt (t-),p Nt (t-)), {a new ,P new ) by increasing characteris- 
tic ratios and call the ordered traits x±, . . . , XN t+ i. Applying Proposition 11.21 
with predators Xi, ...,X]y t+ i, we obtain a globally attracting equilibrium a with 
o~i > if and only if i < m for some m < N t + 1. We then set (ajj(t), = x i 

for all i < m and set (atj(t), Pi(t)) = (0, 0) for all i > m. 

Simulations suggest that we see a growing cloud of coexisting predators 
with some limiting shape and all predators have consumption rates a going 
off to infinity and log(£) going to — oo (see Figures [1] and [2]) . We were not able 
to analyze the two dimensional system so we will specialize to the two cases 
where only a or 5 varies and the other remains fixed. 

1.3. The Alpha Predator Evolution Process 

In this section, we assume that 5 = 1 remains fixed and allow for mutations 
in a. Note that the characteristic ratio of a predator is now 1/a so that 
ordering predators by increasing characteristic ratios is equivalent to ordering 
predators by decreasing consumption rates. We shall prove our results for a 
discrete time version of the Predator EP in which the n th mutation occurs 
at time n. Once this is done, it is straightforward to generalize the result to 
continuous time. 

To more precisely describe the version we study, suppose that at time 0, 
we have a single predator a(0) which can coexist with the prey. If at time n, 
we have N n predators oti(n), ajsf n (n) in decreasing order that satisfy (jl.3p 
with k = N n , then at time n + 1 we choose one of the dj(n), j < N n at 
random, introduce a mutant with trait a new = cij(n) +ell, U ~ Uniform [-1,1], 
and then use ( ll.3p to decide on the state of the process at time n + 1 in the 
same manner as we did for the Predator EP. We shall refer to this process as 
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the Alpha Predator Evolution Process (APEP) and use the following notation 
throughout the remainder of the section: 

• N n = number of coexisting predators at time n. 

• aj(n) = j th largest a amongst all coexisting predators at time n for 
j < N n 

• a mi n(n) = a Nn (n) 

• aj (ri) = a min (n) if j > N n . 

• dj(n) = oij{n) — a min (n) = differences between predator consumption 



• A n = (d 1 (n),d 2 (n), . . .) 
Setting atj(n) = a m i n (n) for j > N n is for convenience so that dj = for all 



Substituting 6 = 1, the condition (11. 3p for coexistence of ax, ajv simpli- 
fies to 



Since a^ja^ > 1, this implies that all the differences aj — o:n must be < r 
so we define S := [0,r] N and let || ■ \\tv denote the total variation norm on 
Aii(S) = space of probability measure on S. We denote by P a the law of the 
APEP started from an initial predator with trait a. 

Theorem 2. iV„ is tight and a m i n (n) — > oo a.s. as n — >■ oo. In addition, there 
exists a measure n e on S and constant a £ > so that 



and a m in{n)/n — > a £ > as n — > oo. 

This result is proved in Section HI The reason for the difference from 
Theorem [1] is that mutant types with traits similar to the resident type can 
always invade because predators only suffer density dependent killings from 
their own type. In the case of prey evolution, this is not the case since inter-and 
intra- specific competition equally affect the prey. 

The key to the proof of Theorem [2] is the observation that as ajv -> oo, 
the condition (jl.4p becomes 



rates. 



j > N n . 




(1.4) 



P° (0) (A n G •) 



Ks(-)\\tv ^0 



N 
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and we can show that the differences A n are asymptotically a positive recurrent 
Harris chain with stationary distribution tt £ . A coupling argument shows that 
the nonhomogeneous chain also converges to tt £ . The linear growth of a m i n 
then follows from a standard result on functionals of positive recurrent Markov 
chains. 

Figure [3] illustrates the tightness of N n and linear growth of a min . Figure 
H] suggests that as the size of the perturbation e — > 0, the spacings between 
traits is 0(e), and the number of coexisting types is 0(l/e). We believe that 
if one converts the re-scaled spacings A n /e into a measure by assigning each 
one mass e then as e — > 0, the distribution of this measure under 7r e converges 
to a deterministic limit in which the density of particles is roughly, but not 
exactly, exponential, see Figure |5j 

1.4- The Delta Predator Evolution Process 

The final evolution process we consider is the Delta Predator Evolution 
Process (DPEP). The DPEP is defined in continuous time and follows the 
same rules as the Predator EP except that all predators have fixed a = 1 and 
we only allow for mutations in 5. For convenience, we assume that e = 1 and 
define Xj(t) = — log Sj(t). We also set N t = the number of coexisting predators 
at time t. Note that since Si(t) < 5 2 (t) < ■■■ < S^f t (t) by definition of the 
DPEP, we always have X x (t) > X 2 (t) > ■ ■ ■ > X Nt {t) so that X max (t) = X x (t) 
and X min (t) = Xjsr t (t) give the traits of the most and least fit predators, 
respectively. Furthermore, (11 .3p implies that 



for all t > 0. 

To get started in the analysis of the DPEP, our first step in Section is to 
prove a simple result which already shows that the behavior is much different 
from the APEP. 

Lemma 1.1. As t — > oo, N t — > oo a.s. 

Our next result describes the asymptotic rates at which the smallest and 
largest predator traits increase. In what follows, we let St be a random walk 
starting at that takes jumps at rate 1 uniform on [—1, 1]. The theory of large 
deviations tells us that 



exists and can be calculated in terms of the moment generating function of St- 




(1.5) 



A(x) = lim - log P(S t > xt) 
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Theorem 3. X max (t)/t — > a and X miri (t)/t — > b a.s. as t — >■ oo where a ~ 
0.9053 an<i 6 ~ 0.5667 satisfy the equations 

A(a) = -1, A(6) = -l + 6. (1.6) 

Furthermore, we have liminf (1/t) log N t > b a.s. 

We will prove Theorem |3] in Section [5j To explain why it is true, it is con- 
venient to adopt the perspective that the different predator types correspond 
to particles and their traits correspond to positions on the real line. Let Z t be 
the branching random walk in which particles give birth at rate 1 and their 
offspring are displaced by an amount uniform on [—1, 1]. A result of Biggins 
[sj implies that r i; the position of right- most particle at time t in the branching 
random walk, has r t /t — > a and 

1 

-\ogZ t ([xt,oo)) 1 + A(x) 

for < x < a so (1/t) logZ t ([bt, oo)) — >■ b. Since we can construct Z t in such 
a way that all particles in X(t) are in Z t , we must have limsupX max (t)/t < a 
a.s. The definition of b and an argument by contradiction using (11.51) gives the 
upper bound limsup^^ X min (t)/t < b for the speed of the left-most particle. 

To bound lim inf t ^oo X mSuX (t) /t, we consider the following branching-selection 
particle system: at any time t, we have M particles with positions Y^(t) > 
■ ■ ■ > Y^f(t), all giving birth at rate one. Whenever a new particle is born, 
we reorder and delete the leftmost particle. Using techniques from the study 
of the APEP in Section HJ we could show that Y^{t)/t — > a M , but instead we 
complete the proof of the first result by showing 

lim liminf Y^(t)/t = a. (1.7) 

M— >oo t— >oo 

Nina Gantert has pointed out to us that Berard and Gouere (2008) have re- 
cently proved 

a - a M ~ C(log(M))- 2 

for a related discrete time model in which all M particle split into two and 
only the right-most M are kept. This confirms a slow rate of convergence, 
which was predicted much earlier by Brunet and Derrida [8j, and which we 
observed in our numerical attempts to verify the limit in (11. 7p . see Figure [61 

To bound liminf^oo X min (/f:)/t, we study the branching random walk with 
killing at — K + jt. Our result given in Lemma 15.41 is a cousin of a re- 
sult of Kesten j25j for branching Brownian motion on [0, oo) where during 
its lifetime, each particle moves according to Brownian Motion with drift 
fi < and variance a 2 , all particles die at rate c and give birth to a mean 
m number of offspring upon death with particles killed when they hit 0. 
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Kesten's result states that the system has positive probability of survival when 
/i < fi = (2cr 2 c(m — l)) 1 / 2 (Theorem 1.1, (1.6)), and in this supercritical case, 
if we start with one particle at x, then for every interval J, Z t (I)/E x Z t (I) — > W 
a.s for some finite random variables W (Theorem 1.1 (1.5)). However, Kesten's 
efforts are concentrated on the exotic behavior in the critical case fi = /io, and 
he says "so far we have only an ugly and complicated proof of the growth re- 
sults in the supercritical case, and we shall therefore not prove Theorem 1.1." 
In section 5, we show that using ideas of Biggins [i| it is easy to prove results 
for (1 /t) log Z t ([ct,oa)). 

The result limsup^^ X min (t)/t < b implies that if T is large and we start 
the branching random walk with one particle at X max (T) at time T then all 
of the particles in the branching random walk with killing at (b + e)t are 
present in the Xi(t). If X min (t) is too far to the left then we would contradict 
fll.5p . The last part of the proof suggests that most particles are near X min (t). 
Simulations (see Figure [7]) further suggest that: 

Conjecture. If we put mass exp(— X min (t)) at Xi(t) — X min (i) then this mea- 
sure converges to a deterministic limit, which again is roughly but not exactly 
exponential. 

However, pro ving; this seems to be a difficult problem. Recently, Durrett and 
Remenik [l7|] have proved convergence of the toy model to the solution of a 
free boundary problem as M — > oo. 

The final conclusion liminf t _ >00 (log N t )/t > b follows from the result for 
X m i n (t) and the proof of Lemma 11.11 Since the result comes from replacing 
ffL5]) by e ~ XN ^ ((3 + N) < r, it seems unlikely that b is the right constant, but 
finding the right constant would require proving the conjecture. 

The proof of liminf 4 _ 5 . 00 X max (t)/t > a also leads to the following result 
regarding the limiting behavior of the most fit predator in the APEP as the 
mutation radius e — > 0. 

Corollary 1. If we run the APEP in continuous time and let 

a* = lim a e m3X {t)/t, 

t— >oo 

then \im £ ^oa* = a. 

2. Prey Evolution 

In this section we will prove Theorem [TJ We first establish Proposition 11.11 
as a consequence of Lemmas 12.11 and 12.21 below. Note that these results only 
cover M < 3, but the proof of Lemma 12.21 can easily be extended to show 
that coexistence for M > 4 is not possible. Our results rely heavily on the 
notion of invadability (see Durrett [15|) and results on Lotka-Volterra systems 
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(see Chapter 13 and 15 of Hofbauer and Sigmund |23j) which we shall quote 
as needed. We shall also make use of the notation introduced in Section 11.11 
and assume throughout this section that N = 1 in (II. II) . Therefore, we have 
a single predator, whose density we shall denote by v, with fixed death rate 
5 > 0. 

2.1. Prey ODE Results 

Let i/i = (oti, Pi) G R+. To determine when the predator and a prey with 
trait yi can coexist, we note that if Pi > 1 then in the absence of predators 
the prey reach an equilibrium density 

o?(yi) = (Pi - 1)/Pi. (2.1) 

If the prey are in equilibrium then the predators can increase when v is small 
if 

ai<r?(yi) - 5 > 0. (2.2) 
Using the formula for a®(yi), we see that this holds if and only if ati > 8 and 

Pi > > 1. (2.3) 

a>i — o 

We call this set of (a>i, Pi) the viable region for prey and label it V. See Figure 
[8] for an example. 

Algebra shows that when (12. 3p occurs, there is a predator-prey equilibrium 
u l (yi) with coordinates 

i, x (#L-l) + ai5 Xl N (ft - l)ai - Pi5 

Pi + <x\ pi + 

A second prey type with trait y 2 = («2, P2) can invade the first prey and the 
predator in equilibrium when 

P 2 (l-a 1 1 (yi))-l-a 2 a 1 2 (yi)>0. 

By interchanging the roles of yi and y 2 we get the condition for the first prey 
to invade the second prey and predator in equilibrium. If both prey traits are 
viable and the two invadability conditions hold, then it is shown in Section 



7.1 of [15J that there is coexistence in the ODE, i.e., the three densities stay 
bounded away from 0. That the densities actually converge to a positive 
equilibrium in this case is the result of Lemma [2 . 1 1 b elow . 

Following [l5[ , we use >~ for "invades" (prey j can invade prey 1 , . . . , j — 1 in 
equilibrium if its density will increase whenever 1, j — 1 are in equilibrium 
with the predator and a small initial density of fs is introduced). Using the 
new notation and defining 

F(yi,y 2 ) = P 2 {1 - o\(yi)) - 1 - a 2 a\(yi), 
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we have 2 >- 1 if and only if 

y 2 G {y : F(y 1 , y) > 0} =: C yi 

and 1 >- 2 if and only if 

y 2 G {y : F(y, yi ) > 0} =: U yi . 

We call the boundary curves L yi = {y : F(yi,y) = 0} and Z7 yi = {y : 
F(y,yx) = 0} the invadability curves and note that we have the formulas 

(a, /3) eU yi <=> /3 = 0(2/1, a) 
(a, (3) e L yi f3 = h(yi, a) 

where 

. N (#l - l)a 2 + (ai - 0i5)a + ft 

9(Vl,<X) = 7— 7 T\ 

1 + «i(a — 0) 

and 

= -j u — T 

1 - 0-1(2/1) 

are well defined provided y±, (a, (3) G V. Calculus shows that the curve U yi is 
tangent to the curve L yi at y\ and we let N(yi) denote the corresponding unit 
normal vector: 

M(y 1 )=c(-a l 2 (y 1 ),l-a{(y 1 )) (2.4) 

where c is chosen to make the length one. The situation is depicted in Figure 
[HJ Lemma 12.11 describes the set of all possible ecological outcomes based on 
this splitting of type space. 

Lemma 2.1. Let y^ y 2 G V with (3i ^ (3 2 and suppose that (ui(0),u 2 (0),v(0)) G 
T^i- Then one of the following must be true 

(a) y 2 G C yi R U Vl in which case the solution to ( 11. 1 It converges to a unique 

positive equilibrium a 2 {y u y 2 ) = (of (y h y 2 ), of(j/i, y 2 ), ^KVh V*)) e r Ji- 

(b) y 2 G £ yi , 2/2 ^ W w in which case the solution to ( 11 . lit converges to the 

equilibrium (0 , a\(y 2 ) , al(y 2 )) . 

(c) y 2 £ Wj/i; 2/2 ^ £ yi in which case the solution to ( II. ip converges to the 

equilibrium {a\{yi), 0, c^j/i)). 

Since inserting a mutant with the same birth rate as the resident is zero, 
the condition ft 7^ ft does not impose any additional restrictions and saves us 
the headache of dealing with a scenario in which we have an infinite number 
of equilibria. 
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Before beginning the proof, we note that setting 113 = v , we can rewrite 
(II. ip in Lotka-Volterra form as 

— = UiiTi + (Au)i) 

where fj = (3i — 1, i = 1, 2, r 3 = —5, and 



.4 



-A -Pi -at 
-f3 2 -(3 2 -OLi 

OL\ «2 "I 



If fi(u) = dui/dt, then the components of the Jacobian are given by 

dfiiu) 



duj 



5ij(ri + (Au)i) +u i A ij . (2.5) 



Let J p denote the value of the Jacobian matrix J = (dfi/duj) evaluated at p. 
Following [23[ (see pages 155, 159), we shall say that an equilibrium point p 
for ( II. ip is regular if det( J p ) 7^ and saturated if + (Ap)i < for all i. Note 
that if we have an equilibrium p 6 T^, then p is trivially saturated since in 
this case, r*j + {Ap)i = for all i. More generally, an equilibrium point a with 
<7j = for all i e J is saturated if it cannot be invaded by types i e /. In 
particular, (o - ^?/!), 0, cr^T/i)) and (0, a{(?/ 2 ), cr^^)) are saturated exactly when 
?/2 ^ £ yi and y 2 ^ respectively. The assumptions yi, y 2 G V and /3i , /3 2 > 1 
imply that cr°(?/i), a (1/2), and the origin are never saturated. It is easy to see 
that since {3± ^ 02, there can be no other possible equilibria a ^ HlV 
Let 

ind{p) = sign(det(J p )) 
denote the index of a regular equilibrium p. The index theorem for Lotka- 



Volterra equations (13.4.4 in 23]) tells us that if all saturated equilibria p are 
regular, we must always have 

ind(p) = (-1) 3 = -1. (2.6) 

p:p saturated 

The key to the proof will be showing that (11.11) has a unique saturated fixed 
point in all three cases (a) - (c). However, since it is not always true that a 
unique saturated fixed point is globally attracting (see pg. 195 in [23]), we 
need to work a little bit harder to get the result. To ease notation, we shall 
let Fi denote the face in r 21 on which Ui = and E^j denote the edge where 
Ui = Uj = 0. 

Proof. Suppose first that we are in case (a) so that (^i(?/i),_0, 02(2/1)) an d 
(0, o~l(y 2 ), o-l(y2)) are not saturated. Then Theorem 7.1 in |l5| implies we 
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have coexistence and hence by Theorems 13.3.1 and 13.5.2 in [23], (11. ip has a 
unique regular equilibrium a 2 G T^i- To show that it is globally attracting, we 



will show that all eigenvalues of J a 2 have negative real parts. The conclusion 
that a 2 is globally attracting on rJi then follows by Theorem 15.3.1 in [23f 
(A is Volterra-Lyapunov stable with di = a 2 ). The Routh-Hurwitz (R-H 
conditions (see pages 702-703 of Murray |27j for the version used here or 
for an elementary proof) tell us that all eigenvalues of J a 2 will have negative 
real parts if (1) trace(J CT 2) < 0, (2) det(J CT 2) < 0, and (3) 

det(J cr 2) > trace (Jo-2)E 2 

where £2 is the sum of the 2x2 principal minors of J a i. But since r + Aa 2 = 0, 
substituting a 2 into (12.51) yields 



J a 2 



-a 2 (3 2 



a 2 a 2 



-a\oti 

-(j\ot 2 

-a 2 



so that the first R-H condition is obvious and the third follows from a simple 
algebraic calculation. The second condition follows from (12.61) since a 2 is the 
unique saturated equilibrium point for (II. ip and is regular. 

Suppose now we are in case (b) so that (0, 01(2/2), 02(2/2)) is saturated, but 
{p\ (yi), 0, a\ (|/i)) is not saturated. If we let = (0,01,02) = (0, 01(2/2), 02(2/2)), 
then 

" n + (Ad)i 

J a = -Ol/?2 -0l/?2 -Ol«2 

2 «1 02«2 —027 

The assumptions in (b) imply that 2 >- 1,3 and 1 )f 2,3 so we must have 
T\ + (Aa)i < 0. Therefore, 

det(J CT ) = (n + (Ar)i)(/3 2 7 + aDa^ < 

implying that a is regular. If coexistence was possible, then as in the proof 
of (a), we would have a regular equilibrium p e T^i and (11.11) would have 
exactly two regular, saturated equilibria, violating (12.61) . Therefore we know 
that Ui(t) — > as t — > 00 for some i — 1, 2, 3. But since we have the invadability 
conditions, 1 >- 0, 3 >~ 1, 3 y 2, 2 >- 1, 3, and 1 )f 2, 3, the proof of Theorem 
7.1 in [l5| tells us that there exists a repelling function for the set 



F = F3 U Ei 2 U F 2 

and therefore we know that (II. ip must leave through F±\F on which a 
is globally attracting (Lemma 5.0 in [HI]). The proof of (c) is identical after 
interchanging the roles of y\ and y 2 . □ 
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Our next result describes the possible outcomes of adding a new species 
when two are already coexisting. 

Lemma 2.2. Let 2/1,2/2, 2/3 EV, y 2 E C yi nU yi , fa ^ f3 2 , and 

(tii(0) ) u 2 (0) ) « 3 (0) l «(0))6lt 1 . 
Then one of the following must be true 

(a) 2/3 G C yi fl C V2 and 2/3 ^ W yi U W<, 2 in which case the solution to (11.11) 

converges to (0, 0, a{(y 3 ), a\(y z )). 

(b) 2/3 E U yi fl W y2 and 2/3 ^ A/i U £j/ 2 ^ n which case the solution to (11. ip 

converges to (07(2/1, 2/2), of (2/1, 2/2), 0, af(yi, y 2 )). 

(c) Neither (a) nor (b) is satisfied in which case either Ui(t) — >■ or u 2 (t) — >■ 

0. 

Proof. It is easy to check that if we have M = 3 in (II. ip . relabel the predator's 
density as M4 and write (ll.ip in Lotka-Volterra form, then det(A) = so that 
there can be no coexistence of the three prey types by Theorem 13.5.2 in [231 ] . 
Therefore, at least one of the types dies out. Which one(s) can be sorted out 
using the same idea as the proof of Lemma 12.11 □ 

To see why we do not need to be concerned with case (c), we note that 
because of the tangency of U yi and L yi , the chance of inserting a mutant which 
can coexist with a current resident prey is always of order 



re 

e~ 2 / x 2 dx = 0(e). 



But then the probability of inserting a third prey which can coexist with two 
resident prey is also 0(e) and therefore, we must wait 0(l/e 2 ) time steps until 
the first time we encounter the situation in Lemma I2T21 (c). Of course, we will 
still get convergence to an equilibrium that can be determined as in Lemma 
O provided & ^/3i,/3 2 . 

2.2. Proof of Theorem^ 

Throughout thus section, we shall use C = Ct to denote a positive constant 
which depends on T and may change from line to line. Write 

where /(2/1, ®) = 0(2/1, a) — h(yi, a) > 0. We also define 



f(t) 



F £ (Y((t)) if lCf(t) = 

max(F £ (Y7(t)), F £ (Kf (t))) if Y£(t) + 0. 
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For Yi(t), Y£(t) G V, p e gives the probability of inserting a new type that 
can coexist when only one resident type is present and an upper bound on 
the probability a new type can coexist with one of the resident types when 
two resident types are present. Note that the tangency of g(yi, ■), h(yx, •) and 
Taylor's theorem imply that 

e- 1 F%y 1 ) -> a p(y 1 ) (2.7) 
for any j/i G V and some constant c > where 

d 2 (g-h) 

W) = g ft 2 (^^l) 

is a continuous function of j/i. Therefore, if Y£(t), Y£(t) G K, K bounded, 
then there exists C# > so that p £ (t) < Ck£- 

Choose an open, bounded set K\ C V and a compact set if 2 C Ki so 
that 2/1 (t) G f° r a h t <T. The existence of i^i is guaranteed since on the 
boundary of the viable region, (3 = a/ (a — S) so that the slope of U yi at y\ is 
implying that N{%)\) points straight up, and it is impossible for yi(t) to leave 
the viable region. Let p > be small enough so that K 2 + p C K\ and define 

r = inf {* : Y £ (t) ^ (K 2 + p) x (K 2 + p)} 

as the first time Yf(t) or F 2 £ (t) leaves K 2 + p. Then p e (t At) < Ce, Vt < T. 
Since mutations occur at rate 1/e, it follows that the expected number of times 
before T At that there is one prey type and an inserted type coexists is < C. If 
two prey types coexist, Lemma l2~2l implies that with probability > 1 — Ce, the 
next time a new type is inserted in C yi fl C y2 , it will replace the two coexisting 
types. Therefore, if e is small, the amount of time during which two types 
coexist is approximately exponential with mean 2e so that 

|0 < T : Y£(t At) > 0}| -»■ (2.8) 

a.s. as £ — > and hence, we can ignore these isolated episodes when studying 
the evolution of Yf(t A t). 

When there is no coexistence, mutations in the direction of L Vl leave the 
resident type unchanged and mutations in the direction of M{y\) replace the 
resident so that the infinitesimal mean of Yf(t A t) is given by 

KVi) = ^(W) ( 2 - 9 ) 

where the 2/3n comes from the fact that if we choose a point at random from 
the upper half of the ball of radius 1 in the (a, 0) plane, then the (3 component 
has density (4/7r)yI — /3 2 and hence mean 

4 f 1 1 4 

- / P^Wdy = — . 

7T J 37T 
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(I2.9P then follows by noting that choices from the half of the ball above L yi 
occur with probability 1/2. It is clear from the scaling that the entries in the 
infinitesimal covariance are of order e and therefore, the infinitesimal mean 
and covariance of Y({- A r) converge to b(yi) and a(yi) = respectively. Since 
b is Lipschitz continuous, the martingale problem for (a, b) is well posed so 
that convergence of Yf(- A r) to y\ follows from Theorem 7.4.1 in Ethier and 



Kurtz [18j]. But then we can choose p small enough so that P(r < T) — > 
and we obtain (11. 2p . 

It remains to prove that 

N* = \{8<t:Y 2 e (a-)=0, Y 2 %s)^0}\ 

converges to a nonhomogeneous Poisson process. Since F 6 (Y^(t)) gives the 
jump probabilities for Nf when Y 2 £ (t) = 0, the compensator for Nf is given by 

A\= \\ {Yi{s)=Q} s- x F\Yt{s))ds. 
Jo 

(EZD, (ESD, and (OD then imply that 

A\ ->■ m(t) = [ ccj)(yi(s))ds. (2.10) 
Jo 

m(t) is continuous and deterministic so we conclude from Theorem 1 in Brown 
[3] that iV £ N where N is a nonhomogeneous Poisson Process with mean 
function m(t). 

□ 



3. Multiple Predator ODE facts 

The goal of this section is the derivation of Lemmas 13.11 and 13.21 which to- 
gether imply Proposition 11.21 The first result gives the algebraic condition for 
existence of positive equilibrium densities and the second proves convergence 
to equilibrium. See Section 11.21 for relevant notation. 

In the absence of predators, the prey have equilibrium density a = r/(3 
where r = (3 — l>0by assumption. Suppose we wish to find a positive 
equilibrium a k = (a*, of , • • • on the face 

ri,fc = {v e r 1>JV : v k+1 = ■ ■ ■ v N = 0}. 

Then, solving the equations atj — 6j — a*- = for of, j = 1, k we obtain 

a* = a k J (x 1 ,... , x k ) = ajdQ - 6j 
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and substituting these expressions into the equation r — (ia^ — EjLi a j al j = 
yields 

k k 
3=1 3=1 

We can conclude that 

for all /c > 0. To determine when <7* > 0, 1 < j < k, write S k = Yli=i a l an d 

aj(T ° = J+s h • 

Adding 8j — 5j({3 + S k )/(f3 + S k ), the above 

af - 06 j + )J ( .. ; (n,nA - afSj) 



Si + 



5j + 



(3 + S k 

k-l 



(0 + Y.i +j at?)(atjcr (x h ...,x,-u -,x k ) ~ Sj) 



since (/3 + a i) a o~ 1 = r + Ei^y a iSi- From this it follows that o k will be 
positive if and only if 

o-g^^i, a^_i, x j+1 , x k ) > ij. (3.1) 
where £j = 8j/cxj is the characteristic ratio for predator Xj. 

Lemma 3.1. Suppose xi, x k are ordered by increasing characteristic ratios. 
Then a k > if and only if 

k 

/34 + J>?(4-^) <r (3.2) 

3=1 

and if (13.21) is satisfied, then a, > for all j < k. 
Proof. (13.11) implies that o\ > if and only if 



k <<T*- 1 (x h ...,x k -i)= n JZk\ i3 , (3-3) 



where we have used the definition of £ 3 - = <5j/<X/ on the right. Multiplying both 
sides by the denominator of the right and then rearranging terms, we obtain 
(13. 2 p (since £ k — £ k = 0). Now suppose that (13. 2p holds. Then since £ k > £j, 
for all j = 1, k — 1, (13. 2p also holds if we replace £ k by £j, j < k on the 
left and reversing the algebra used to derive (13 .2p from (13.31) shows that this 
is equivalent to > 0, proving the result. □ 
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Lemma 3.2. Suppose we have a collection of predators xi,...,xn ordered by 
increasing I 's and let k < N be the largest integer for which (13 ,2p is sat- 
isfied (with the convention that k = if (13. 2 1) fails for all k < N). Then 
a = (erg , of, . . . , cr£, 0, . . . ,0) is a globally attracting fixed point on T^ N with 
Lyapunov function 



N 

Vi- 

i=l i=k+l 



V(u,vx, ...,v N ) =u- <t{ logw + ^(Ui - o-flogUj) + 

Proof. Differentiating V yields 

dV k n 

j=l i=fc+l 



i=l i=fc+l 

fc jV 2V 



i=l i=ft+l i=fc+l 

If (w, Ui, Vn) = o~, this expression is and otherwise it is < since Lemma 
Oand (JO) imply that 

A > 4+1 > 

for all i > k + 1 so that all terms on the left are negative. 

□ 

4. Proof of Theorem [2] 

In this section, we prove Theorem [2] for the APEP as defined in Section 
11.31 and use the notation defined there. We also define the Markov chain 
Y n = (a min (n), A n ). 

If 5j = 1 for all j, (13 .2D with k = N can be rewritten as: 

N R 

Y] — («i - a N ) <r - —. (4.1) 

Our first step is to show 

Lemma 4.1. T/ie sequence N n is tight. 

Proof. Define the sets A m = [0, r] m x {0} N , for m G N. Then A n G A m if and 
only if N n < m. Let M = M(r,e) = [~-^] be the smallest integer > 4r/e and 
suppose that F„ = t/e M + x 5. From (14. ip . at most M of the aj(n)'s can be 
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> a m i n (n) +e/4. With probability at least 1/4 M , the next M mutants will be 
inserted to the right of a m i n (^) + e/2. But then none of the predators to the 
left a m j n (n) + e/4 can be in the coexisting set at time n + M because otherwise, 
by the definition of the APEP, any predator with a > a m \ n {n) + e/2 would 
also be in the set, and since there are at least M such predators, 

oo 

dj{n + M)> M (e/2 - e/4) > r 

i=i 

contradicting (14. ip . Therefore, we have the uniform lower bound 

P(A n+M eA 2M \Y n = y)>4- M (4.2) 

which holds for all y e M + x S. Since this bound is uniform in y, tightness 
follows. □ 

Lemma 4.2. As n — > oo, the marginal transition probabilities for A n : 

p a {A,-) :=P(A n+1 e -\Y n = {a,A)) 

converge in total variation to the transition probabilities for a time homoge- 
neous Markov Chain X n with transition probabilities p(A, ■). 

Proof. Suppose that Y n = (a, A). Since < ay(n) —a m i n (n) < r for all n > 1, 
we can see that as a — > oo, (14. ID simplifies to 

N 
3=1 

This implies that, in the limit, the differences evolve according to the following 
algorithm: pick a species 1 < k < N n at random, insert a random mutation 
in (dk(n) — e,dk{n) + e), and then modify the algorithm in Theorem 11.21 to 
use (14. 3[) instead of (14. ip with the rule that we shift the differences before 
calculating the sum if the new insertion is left of 0. □ 

Our next result concerns the limiting behavior of X n . Writing x instead of 
A for the vector of differences, we set 

p{x,A) = P{X n+1 eA\X n = x). 

Lemma 4.3. X n is a positive recurrent, Harris Chain and hence, has a unique 
stationary distribution ir. 

Proof. Following the arguments in Athreya and Ney [2j, it suffices to show 
that there exists a "regenerative" set A C S satisfying: 



21 



(CI) P x (ta < oo) = 1 for all x £ S where ta is hitting time of A. 



(C2) There exists a probability measure p on A, A > 0, and k £ N so that 
p K (x, B) > Xp(B) for all x £ A, B C A. 

The same calculation that led to (14. 2D shows that y4 2A/ satisfies the con- 
dition in (CI), but (C2) may not hold for this set. We therefore define a set 
G (for good) that will be reached from A 2M with probability 1 and satisfies 
(C2). To this end, let 

1 + sup { k : y * j = ^ J" - < 2r/e 




and choose rj small enough so that 

k 

^j(e/2 + 77)<r. (4.4) 

3=1 

Let G = {di — di+i £ (s/2, e/2 + rf) for i < k and di = for i > k}. In other 
words, d £ G corresponds to k types coexisting with a's that have spacings 
between e/2 and e/2 + rj units apart. 

The first step in showing that (CI) and (C2) hold for A = G is to show 
that if Xq = x £ A 2M , then we can get to A in k steps by the following 
path: first, we choose di (the predator with the largest values of a) as our 
mutating predator at time 1 (which happens with probability at least (2M) _1 ) 
and then choose a mutant type g± in (d\ + e/2, d\ + e/2 + rj) (which happens 
with probability rj/(2e)). The next time step, we choose g± as our mutating 
type (which happens with probability at least (2M + 1) _1 ) and then mutate 
to g>2 £ (gi + z/2, gi + e/2 + rj). If we continue for k steps, then each gj, 
1 < j < K will be at least as big as di + je/2 so that by (14. 3p . no member 
of the coexisting set at time will remain at time k. Furthermore, by (14. 4p . 
the shifted set d'j = g K -j+i — gi, 1 < j < k will satisfy (14. 3 \ and therefore, 
X K £ G. It is clear from the construction that we have 



p K (x,G)>[^-^ (4.5) 

F v ' ; ~ \2e{2M + K) J K J 

To prove (C2) holds, we first consider cylinder sets of the form B = {di — 
di + i £ Bi C (e/2, e/2 + 77) for i < k and di = for i > k}. Then if x £ G, 
taking the same path that led to (14. 5 j) yields the lower bound 

p K (x, 5) > C _J__] K . (4.6) 
y V 7 - (2e) K - 1 V 2M + K / 
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If we let p = Lebesgue measure on G normalized to be a probability and 
recall that the Radon-Nikodym derivative dp K (x, -)/dp(-) evaluated at a general 
measurable set B can be written as the limit of p K (x, Bk)) / p(Bk) where Bk is 
a sequence of cylinder sets, (C2) follows. 

To check positive recurrence, we let ta be the first hitting time of our 
regenerative set G. (14. 2 p and ( 14. 5 p tell us that there is a positive constant 
9 = 9(r, e) so that 

p 2M+K (x,G) >9>0 

for any x £ S. Therefore, we have E x (ja) < (2M + k)/9 < oo, completing the 
proof. □ 

The construction in the previous lemma also yields: 
Lemma 4.4. a min (n) -> oo a.s. as n —¥ oo. 

Proof. We can modify the construction in the previous Lemma to show that 
there exist constants K, J > 1, p > so that 

P{ ai {{n + l)K) - a x (nK) > Je/2\Y nK = y) > p 

for any y £ R + x S and n > 0. Therefore, ai(n) — > oo a.s. by the Borel- 
Cantelli Lemma and the result follows since ai(n) — a m i n (n) < r. □ 

Theorem 4. As n — > oo, ||P a (A n £ ■) — 7r(-)||yy — > for any initial a £ M + . 

Proof. It suffices to prove the result for the subsequences n = vtlk + j for 
< j < k, but then by using the Markov property at time j, it is enough 
to prove the result for n = ran and a general initial distribution. To prepare 
for the proof, recall that one can modify the state space of a Harris recurrent 
Markov chain to have a point ( that corresponds to being distributed on the 
set A according to p with the exact position being independent of the past. 

To prove the result, we will construct a process (X n , A n ) on S x S so that 
the marginal law of A n is the law of A raK , the marginal distribution of X n is 
7r for all n, and P(X n ^ A n ) — > as n — > oo. Let Lq, U 2 , ■ ■ ■ and Vi, V 2 , . . . 
be independent and uniform on [0, 1]. To begin, let q a (x, ■) = P(A K £ -|A = 
(a,x)) and 

q(x, •) = lim q a (x, •) = p K (x, •) 
by Lemma 14.21 Define the function J n : S x [0, 1] — > S by 

P{J n {x,U n ) £ B) = q amSa ( nK )(x,B). 

Since q auiin (nK)(x, ■) £ Ai\(S) and S is a separable metric space, defining J n is 
possible by Theorem 3.2 in Billingsley [5|. Suppose that X n has distribution 
7T , define Z n+1 = J n (X n , C/ n ) and 



/i n (A) = P(Z n £ v4|a min (nK)) = / q arn . n(nK )(x,A)n(dx) 
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and let (X n+ i, Z n+ i) be a maximal coupling of (X n , Z n ) so that 

P(X n+ i ^ Z n+1 ) = \\/j, n - 7l\\ TV 



(see, for example Thorisson 
we have 



Then from the definition of \i n and (X n+1 , Z n+ i) 



Vn+l = P(^n+1 7^ Zn+l) 

Qa miu (nK)(x, -)^(dx) - I q(x,-)ir(dx) 



->■ 

TV 

as n — > oo by Lemma 14.21 and 14.41 

When {A n = X n }, we set A n+1 = J n {X n , U n ) = Z n+l so that 

X n = A 

n) S Vn+1 

On {X n 7^ A n }, we take A n+ i = J n (X n , V n ). (14.61) implies that g(x, £) > A, so 
it follows from Lemma I4T21 that if a min (nfi) > a then q(a min (nn), x, () > A/2, 
and we have 

P(X n+ i = A n+1 |X„ 7^ A n ) > A/2 
so that if ( n = P(X n 7^ A n ), then 

Cn+i < (1 - A/2)C„ + 77 n+ i. 

Iterating, yields the inequality 



n+1 



Cn+1 < J> - V 2 ) 



n+l— i 



(4.7) 



i=i 



Since |1 — A/2 j < 1 and r]„ — >■ 0, the right hand side of (14. 7p must also go to 
zero which yields 

\\P(A nK G •) - vr(-)||rv < P{X n ^ A n ) = r? n -> 0. 
completing the proof. □ 

It remains to prove the result on the linear growth of a m i n (n). Since atj — 
®min < r , it suffices to establish this for a max . To do this, we look at the chains 
Z n = (X n , U n , V n ) with U n uniform on [0, 1] giving the index k = \N n U n ~\ of the 
value to be mutated, and V n independent uniform on [— e, e] giving the change 
in the value due to mutation. It is clear that the distribution of Z n will converge 
in distribution to the product measure n = it x uniform[0, 1] x uniformf— e, e] 
so that if we let f(Z n ) = a m3uX (n) — a max (^ — 1) be the amount shifted at the 
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n th step, then / is non-negative and bounded above by e so the strong law for 
functionals of Markov chains implies 



ttmaxO) - a m ax(0) 



n 

m=l 



i - r 

-J2f(Z m )^ J f(x)n(dx) = a (U 



Since / > with positive probability, a > 0. To extend this result to the 
real chain, let (X n , A„) be the coupled chain from the proof of Theorem H] and 
define D n = 1 if X n ^ A n and D n = otherwise. From the proof of Theorem 
HI we can dominate D n by a Markov Chain B n that has 

P(B n+ i = l\B n = 0) = Vn+1 
P(B n+1 = 0\B n = 1) = - 

i.e., we can define the two processes on the same space so that B n > D n for 
all n. Coupling B n with a homogeneous chain B£ that has P(B^ +l = 1\B? = 
0) = p, P(B^ +1 = 0\B? = 1) = A/2, and stationary distribution n p with 
7r p (l) — p/(p+ A/2) and recalling that r\ n — >■ 0, it follows that 

limsup - V" -D m < lim sup - V" 5 m < — -^-777- 
n^oo n ^— ' n^oo n ^— ' p + A/2 

m=l m=l ' 

Since this holds for any p > 0, we must have 

1 n 

lim — \ D m = 

n— >oo 77, ^— — ' 



n 

m=l 



and the desired result now follows from (14.81) and the fact that < a m3uX (n) — 
a mffl (?i — 1) < 1 for all n > 0. 

5. Proof of Theorem [3] 

In this section, we prove Theorem [2] for the DPEP as defined in Section [L3l 
and use the notation defined there. Since one of the keys to deriving our results 
will be comparison with a branching random walk, we continue adopting the 
perspective that Xj (t) refers to the position of particle j on positive half line. 
Note that if we set aij = 1, k = N in ( 13.21) . we obtain the condition for 
coexistence: 



>N 



(»s(-a) 



<r. (5.1) 



Proof of LemmaUUl Let X^t) = -log(^(t)) and X^t) > ■ ■ > X M (t) be the 
rightmost M particles at this time. It should be clear from (15. ip that if 



r XM{T \(3 + M) < r (5.2) 
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then we will have N t > M for t > T. Let y = -log(r/(/3 + M)). The right 
most particle is increasing in t. Since the number of particles changes by < 1 
each time and Ylm=i V m = 00 the right-most particle gives birth to the right 
of its current position plus 1/2 infinitely many times. Thus at some time T, 
we will have a point > y + M. Since \Xi(t) — X i+ i(t)\ < 1 and points are only 
erased when (15. ip fails, (15. 2p follows. □ 

5. 1 . Asymptotics for X max 

For the remainder of the paper, we let Z t be a branching random walk 
started from one particle at 0, in which particles give birth at rate 1 and 
displacements are uniform on [—1,1]. It is well known that the mean measure 

EZ t {A) = e'PiSt e A) (5.3) 

where S t is a continuous time random walk that jumps at rate one and takes 
step uniform on [—1, 1]. If we let (p(8) = (e e — e~ e )/28 be the moment gener- 
ating function for the displacements, then 

00 4.n 

Ee es * = V e -'-^(0) = exp(t(#0) - 1)) 
^— ' n! 

n=0 

Chebyshev's inequality implies that if 9 > 

P(S t > xt) < exp(-t(6x - (f)(6) + 1)) (5.4) 
and standard large deviations results imply that for x > 0, 

- log P(S t > xt) A(x) = - (sup{9x - <p(9)} + 1 ) (5.5) 

* \6>0 J 

where A(0) = and A is strictly decreasing on [0, 00). 

Biggins 0] , Theorem 2 shows that the right-most particle in the branching 
random walk Z ma _ x (t)/t — >• a a.s. where a, defined in (\l.6\1 . is the smallest 
x > such that A(x) < —1. Since the particles Xi(t) in our evolution model 
are a subset of those in the branching random walk, we have 

limsup Xi(t)/t < a. 

t— >oo 

The remainder of this section is dedicated to the proof of the lower bound 

liminf X x {t)/t > a. (5.6) 

t— >oo 

By Lemma ll.l[ we know there exists some time T so that N t > M for 
t > T. By the proof of Lemma 11.11 we can take T to be the first time 
e -x M (T)^ _|_ < r ^ -^i,^ j s a stopping time, so the future behavior of the 
process is not affected. 
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Lemma 5.1. If we start the toy model at time T with positions equal to the 
rightmost M particles at this time Xi(T) > ■ ■ ■ > Xm{T), then the Xi(t) 
and Y/ 4 (t) can be defined on the same space so that Xi(t) > Y i (t) for all 
l<i<M andt>T. 

Proof. Couple the birth times of X^it) and Yj M (t) and the displacements of 
their offspring. Recall that if a birth from Xk(t) with k > M lands to the right 
of some Xi(t), i < M, we renumber the Xi and put them in decreasing order. 
Births of particles from Xk(t) for k > M may cause the X's to get ahead of 
the K's, but coupled births for i < M cause the vectors of X's and K's to 
move in parallel and the desired comparison follows. □ 

For our next comparison consider the branching random walk started with 
one particle at K x (0). Let Tk be the time of the kth birth, with To = 0, 
and for t G [T A _ 1; T fc ) let tf(t) > C*(*) > • • > C*(*) be the locations of the 
particles present. 



Lemma 5.2. We can couple the branching random walk and the toy model so 



that for t G [T fc _i, T k ), Y M (t) > (-(t) for 1 < j < k and k < M. 



Proof. Couple the birth times of Cj{t) an d Y^ M (t) for j < k and t G T k ], 
i.e., there will be no births in (Tk-i,Tk) and the same particle will give birth 
at time Births of particles from Yj(t) for j > k may cause the K's to get 
ahead of the C's, but coupled births for j < k cause the vectors of £'s and K's 
to move in parallel. □ 

Lemma 5.3. Let Bm be the time of the Mth birth in the branching random 
walk. 

hminf^> EZ ™L {Bm) osM-ioo 
t^oo t EBm 

Proof. Let Tj^i = Tk where Tk are as in Lemma 15.21 and for j > 1, let Tkj, 
k < M denote the time of the k th birth in a BRW started with a single particle 
at Y^' 1 {T M j_i) at time T M j_i and let £mj denote the position of the rightmost 
particle at time Tmj- Repeatedly applying the comparison in Lemma l5\2l and 
letting yields 

Y?® > E j: T M , 3 <t( Y i M ( T M, j+ i) - Ki M (T M ,j)) 
t ~ t 



t 



But the time intervals Tmj+i — Tmj are iid with mean EBm so the first part 
of the result follows from the renewal theorem. To prove the second part, we 
note that Biggins' result implies 

Zul^Bm) I Bm —t- a almost surely. 
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Since Bm = + - ■ - + £m where the & are independent exponentials with mean 
it is easy to see that B M /EB M — > 1, so 

Z ma , x (B M )/ EB M — >• a almost surely. 

Therefore, the result will follow from the dominated convergence theorem if 
we can show that 

By the Cauchy Schwarz inequality, it suffices to show 

E sup — < oo (5.7) 

\t>i t J 



and 



To prove (15. 7ft . we note that (15. 3 h and (15.41) imply that 

P(Z max (t) > xt) < e' (1+A(a;)) 

and since A is concave with A(0) = and A (a) = —1 with a < 1, it follows 
that for x > 1 

P(Z max (t) > xt) < 

Now if Z max (t)/t > 2x for some t, then since Z milx (t) is non- decreasing, we 
must have Z max (s)/s > x for some s£ [t, t + 1] and therefore, summing over 
all integers t from 1 to oo, we see that if x > 2 



P supZ max (0A>2x <e 



which proves (15. 7p . To prove (I5.8p . we note that -E-Bm = Z~2iti V* an d 



Sexp(0B M ) = H — 



for < 9 < 1, so by Chebyshev's inequality, 



9/i 



( M 1 - \ 

P(B M > yEB M ) < exp \-9y - ~ J2 h ^ ~ 9 / i ) 

\ i=i 1 i=i / 
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Taking 6 = 1/2 and choosing c so that log(l— x) > — x — cx 2 when < x < 1/2, 
we have 



P(B M > yEB M ) < exp HT -(1 - y) + — 



< Cexp ( — ^-log(M + l) 




) 
) 



C{M + l) {1 - y)/2 



Therefore if y > 3 



oo 



oo 



2(3-y)/2 



^(M + l)^/ 2 ^ 



(v - 3)/2 



M=2 



which yields (I5.8p . completing the proof. 



□ 



(15. 6p follows from Lemma 15.11 and 15.31 which completes the proof that the 
speed of the rightmost particle is a. We shall complete the proof of Theorem 
[3] in the next section by showing the speed of the leftmost particle is b, but 
first we pause to prove Corollary [TJ 

Proof of Corollary [I] Suppose we choose e small enough so that eM(M — 
l)/2 < r. Using the coupling in Lemma 15.21 we can use the particles 
j < k < M, from the branching random walk started at X max to get a lower 
bound on the right-most k < M particles in the predator evolution with fixed 
5. An induction argument shows that the spacings between the corresponding 
particles in the predator evolution are < e at all times. Since we have assumed 
e ^jLi 1 j < r t ne right- most k < M particles are never killed. The remainder 
of the proof is the same as before. □ 

5.2. Asymptotics for X min 

In order to get the speed of the leftmost particle, we will need the following 
result on a branching random walk with killing which is an adaptation of 
Biggins [if, Theorems 1 and 2, which proves this result without killing. 

Lemma 5.4. Let Z t (7, A) denote the number of particles in A under a branch- 
ing random walk with birth rate one, displacements uniform on [—1,1], killing 
to the left of —K + yt, and started with one particle at 0. Then for any c > 7 
on the set of nonextinction 



where 1(c) = 1 + A(c), and the probability of extinction tends to as K — > 00. 



lim - log Zj(7, [ct, 00)) = 1(c) 



(5.9) 
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Proof. Theorem 2 in Biggins along with (15. 3p and (I5.5P yields (15. 9p in the case 
of no killing and since Z t (j, [ct, 00)) C Z t ([ct, 00)), we get the upper bound 
in (15. 91) . To get the lower bound, we recall that to prove the corresponding 
lower bound for the process without killing, Biggins lets Z^ +1 be the points 
at time (m + l)k that are at least kc units to the right of their ancestor in 
Z^ at time mk. \Z^\ is a branching process with offspring distribution [Z\\ 
so (\Z^\) 1 / m — > E\Z\\ on the nonextinction set. Combining (15. 3p and (15. 5p 
implies (1/k) log E\Z^\ — > 1(c) which yields the desired lower bound. 

To extend this construction to the process with killing, let Z^ +1 be the 
points at time (m + l)k that are at least kc units to the right of their ancestor 
in 2!^ at time mk and are not killed by going to the left of —K + jt of 
mk < t < (m + l)k. By construction, all points in are > cmk and we 
have chosen 7 < c so for large m, the killing has little effect and on the set of 
non-extinction we have 

■i]og|z£|-Hogtf|Z*|. 
m 

Using ( 15.31) and (15.5P again gives the desired lower bound. □ 

With this result in hand, we can complete the 
Proof of X min (t) jt b. When X min (t) increases we must have 

Nte -X min (t) > r _ 

Since the particles in X are a subset of the particles in the branching random 
walk, it follows that if X min (t) > (b + e)t, 

Nte -x min (t) < Z( Q( 6 + £ ) f) oo)) e -(^)* _> 

as t — > 00 since 1(c) < c for all c > b. Therefore, lim sxip X m i n (t) /t ^ b a.s. 

To prove that liminf X min (t)/t > b a.s., let c G (b, a) and e > 0. Choose if 
large enough so that the probability of extinction in the branching random walk 
with killing at — k + bt is less than e for all k > K and then take T large enough 
so that Xi(t) > ct, for alH > T (which is possible since limXi(t)/t = a) and 
so that bT > K. Suppose that X min (t) < (b — p)t for some p > 0. Then by 
comparing with a branching random walk with killing at — X\(T) + bt, we have 

Nt 

F(t) := e - x ^ {t) ^(1 - e-^'W/ x ^W) (5.10) 
3=1 

> e - {b - £)t (l - e- {c - b+£)t )Z t (b, [ct, 00)). 
But on the non-extinction set (which has probability at least 1 — e), we have 

lim - log[e- (fe - £)t (l - e- {c - b+£) t)Z t (-f, [ct, 00))] = 1(c) -b + e^e>0 

t 
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as c 4- b and therefore, we must have X min (t) > (b — e)t eventually or there 
would exist a sequence of points £j — > oo for which F{ti) — > oo, contradicting 
(II. 51) . Therefore, P(liminf X m i n (i)/t < 6) < e and since e is arbitrary, this 
proves the result. 

To conclude that liminf t _ >00 (log N t )/t > b a.s., note that if e > then for 
large times there are at least exp((/(c) — e)t) points of X to the left of ct. 
Picking c close to b and e small gives the desired result. □ 
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Figure 1: The five clusters, from upper left to lower right, are the characteristics of the 
coexisting predators in a sample run of the Predator EP after 10 4 , 1.25 x 10 4 , 1.5 x 10 4 , 1.75 x 
10 4 and 2 x 10 4 mutations have occurred. The consumption rates, a, of all coexisting 
predators are plotted on the x-axis and the corresponding values of log/? = \og(S/a) are 
plotted on the y-axis. Parameters: r = 1, a(0) = 3, (5(0) = .45, e = .01. 




Figure 2: Plot of N n = number of coexisting species in the population after the n mutation 
has occurred in the Predator EP from Figure [1] 




Figure 3: Top: Plot of the number of coexisting predators at time n in the APEP with 
e = .01, r = 1, a(Q) = 3. Bottom: Plot of the change in a m i n (n) for the same simulation. 
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Figure 4: The top panel shows the inverse of the average number of species and the bottom 
panel shows the average maximum distance between a's as a function of e for the APEP. 
Here, we have run one simulation for each value of e = .001, .002, ... , .01 with r = 1, 
a(0) = 3 and then averaged out the results of each simulation over the last 25,000 time 
steps to obtain the values for the plotted points. The solid lines are the corresponding least 
square lines. It appears that the number of coexisting species is 0(1/ e) and the maximum 
distance between coexisting types is 0(e) as e — > 0. 
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Figure 5: Plot of the distribution of predator types for a single run of the APEP at time 
n = 50,000 with e = .001 and r = 1. The solid line connects the points (dj/e,e(N n — j)), 
j = 1, 2, . . . N n — 17626. The dashed line gives an exponential approximation. 




0.2- 
0.1 - 



0.5 1 1.5 2 2.5 3 

log M 

Figure 6: Graph of the speed «m versus logM showing slow convergence to the limit 
a w 0.9053 for the finite branching-selection particle system Y defined in Section H~4l 
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Figure 7: Plot of the distribution of predator types for a single run of the DPEP with 
r = 1 at time t ~ 20.25 (after n = 50,000 insertions). The solid line shows the point 
{Xj(t) - X min (t),(N t - j)e- x """ (t *>), j = l,...,N t = 25467. The dashed line gives an 
exponential approximation. 
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Figure 8: Plot of the invadability curves for = (2,4). For the predator, we set 5 = 1. 

The dashed line shows the boundary of the viable region. 
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